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r~! . Abstract 

■ Recent works emphasized the interest of numerical solution of PDE's with wavelets. 
In their works [3, 4], A. Cohen, W. Dahmen and R. DeVore focussed on the non linear 
approximation aspect of the wavelet approximation of PDE's to prove the relevance 
of such methods. In order to extend these results, we focuss on the convergence of 

■ the iterative algorithm, and we consider different possibilities offered by the wavelet 
■^j- \ theory: the tensorial wavelets and the derivation/integration of wavelet bases. We 

also investigate the use of wavelet packets. We apply these extended results to prove 
in the case of the Shannon wavelets, the convergence of the algorithm introduced in 

■ [8]. This algorithm carries out the Lcray projector with divergence-free wavelets. 

O ' 1 Introduction 



Since the end of the 80's, the mathematical Theory of Wavelets has invented new tools 
for numerical simulations. Thanks to the Fast Wavelet Transform, wavelets provide effi- 
cient algorithms including optimal preconditionners for elliptic operators [12]. Recently, 
Cohen-Dahmen-De Vore's articles [3, 4] demonstrated the optimal complexity of wavelet 
algorithms for the solution of elliptic problems, 
jj] ■ These works enhanced the interest for these methods applied to the solution of partial 



X 



differential equations. Wavelet approach also resulted in non-linear approximations [5], 
and in denoising methods [10]. Wavelet methods offer the possibility to regulate both the 
accuracy in space and the accuracy in frequency. 

In the following, we apply the Shannon decomposition to differential operators in order 
to investigate the convergence of wavelet algorithms to solve partial differential equations. 
In the first part, we recall some operator theory basic elements; we indicate how the 
Shannon wavelets and wavelet packets can be used to part the support of the Fourier 
transform of a function. 

In the following part, we resume the works [3, 4] on the wavelet approximation of dif- 
ferential operators and state the correlated theorem of convergence for constant coefficient 
operators. In the frame of Shannon wavelets, we show that this convergence depends on 
the chosen wavelet decomposition (MRA or tensorial). We extend this study to Shannon 
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wavelet packets. After that, we introduce a result on the derivation of biorthogonal wave- 
lets due to P. G. Lemarie-Rieusset [14] that serves to construct new wavelet approximations 
of differential operators. 

In the last two parts, we present explicit examples that are implicated in the numerical 
solution of the Navier-Stokes equations: the solution of the implicit Laplacian (Id— aA) -1 
and the Leray projector P. 

2 Symbol of an operator 

The convergence of the wavelet methods involves the partition of the spectra of the opera- 
tor. This partition is provided by the wavelet decomposition. Hence we need the notion of 
symbol as introduced by Lars Hdrmander in his book The Analysis of Partial Differencial 
Operators [11]. 

We denote by dj the derivation along the variable Xj, and by Dj the derivation — idj. 
For a = (ai, . . . , ay) € N d , we write D a = D" 1 . . . D^ d . Let u denote a Shwartz function 
of d real variables (i.e. u is C°° and fast decreasing: V7V € N, 3A > / Vx € ffi d , \u(x)\ < 
A/(l + \x\ 2 ) N / 2 ). We denote by < • , • > the scalar product either on vectors either in dual 
spaces. 

Then u stands for the Fourier transform of u, i.e. 

u(0 = [ e- i<x 't>u(x)dx 

We also denote by T the isomorphism of L 2 (R d ) given by T : u *—> (27r)~ d / 2 M. 
The inverse Fourier transform is done by: 

u(x) = (2ir)- d / 2 [ e i<x '^ Fu(£)d£ (2.1) 
When we derivate the relation (2.1), it yields: 

D a u(x) = (2vr)- d / 2 f e <<x ' {> ^ a ^u(^)dC 
J^m d 

Thus derivating u by D a consists in multiplying the Fourier transform of u by £ a = 
^l 1 ■ ■ ■ C° ' ■ ^ ne function £ £ a is called the symbol of D a . More generally speaking, if 
a(£) is a C°° function slowly increasing (i.e. such that 3N € N, A > / V£ £ R d , |a(£)| < 
^4(1 + l^l 2 )^/ 2 ), a(D) defines an operator of symbol a(£) acting on the class of the Schwartz 
functions S by 

a(D)u(x) = (2vr)- d / 2 f e i<x ' 5> a(0-M£R (2.2) 

Let's now consider a differential operator P of order m with variable coefficients a a in S, 
P = X^ui=o o- a (x)D a . Then, instead of using the formula: 

Pu = (2vr)- d / 2 f e i<x ^> T{Pu){i)di 
J^eR d 
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where 

m 

T(Pu)(0 = (2vr) d / 2 ]T Ta a *Z a Tu 
|a|=0 

that is no more a multiplication but an integral operator on J-u, we use the formula (for 
x € R d ): 

Pu{x) = (2^ 2 [ [ V a a (x)r | 



that we write 

Pu(x) = (2vr)- d / 2 [ e^'^pfoO^HKR (2.3) 
introducing the "symbol" p(x, £) of P 



\a\=0 

The formula (2.3) gives us the possibility to define the operators p(x, D) of symbol p(x, £) 
that are not polynomials in £. These operators are called pseudo-differential. The functions 
p must verify regularity and increase properties of polynomial type (see [11]). 

We'll remark that the function ^{Pu) is no more the function p{x^)J-u{£ > ) which 
appears in (2.3) since the latter depends on x. 

The following definition of an elliptic operator is given in [11]: 

Definition 2.1 (elliptic operators) From the symbol p(x,^) = X^™|=o a a( x )£, a we ea> 
tract the principal symbol p m (x , (,) = Yl\a\=m a a(%)£, a - A differential operator P is said to 
be elliptic iff 

V^R d , V£eM d \{(0,...,0)}, Pm ( x ,o^o 



In the following, we'll need differential operators applied to vector functions M. d — > M m . 
We denote by u with bold caracter the (multi-variable) vector function u of real variables 
when it has several components. For u having several components, the symbol p(x, £) is a 
matrix: Vx,£ 6 R d , M(x,£) G C nxm . Let A = (Ay) 

i<i<n, i<j<m be a differential operator 

A : (H s (R d )) m -v (H r {R d )) n , with: 

Aij = ^ ^ aij a (x)D 

a 

Its symbol is M = (my)i<i< n i<j< m with 

mij(x,0 =J2a ij:a (x)£ a 
a 

We apply the operator A componentwise as follows: 

m 

(Au)i = ^2 i m i j(x,D)uj 
i=i 

Therefore, the multidimensional symbol can be handled in much the same way as in 1-D. 
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Remark 2.1 As the operator is applied to real functional spaces, its symbol verifies the 
same relation as the Fourier transform of real functions, that is: 

Mi,j mij(x,-g) = my(x,£)* 

where z* denotes the complex conjugate of z. 

3 Shannon wavelet decomposition 

A good reference for the definition of the Shannon wavelets is Mallat's academic book [15]. 
We first briefly recall the construction of these wavelets. 

The biorthogonal wavelets are based on scale filters m and n that provide the 
low-pass filter and the high-pass filter. As the scale function (f(^) and the wavelet ip(^) 
belong to Vq = span{ip(- — k), k € Z}, there exist two sequences (a^) and (b^) such that: 

¥>(|) = ^2a k ip(x - k) , = J2 bk ^ x ~ fc ) 

We thus get after a Fourier Transform: 

with m(0 = \ £ fceZ a k e~^ , n(£) = ± £ fceZ M"** 
The scale function is inferred from the filter by : 

0(0 = 0(0) jl 

3=1 

A wavelet basis {^jk}j,kez with if)jk(x) = 2^ 2 ip(2^x — k) forms a Riesz basis of L 2 (M). 
Similarly, if we denote by H l the set of distribution functions / such that (1+|£| 2 )*^ 2 / € L 2 , 
{2 tj ipjk}j,kez provides a Riesz basis of the Hilbert space iJ*(R) for — ti < t < t2- Where, 
if (cf), if)) is the dual wavelet basis associated to ((f), if)), t\ € M is the maximal number such 
that ip € H tl (M.) and £2 € R the maximal number such that V> £ i?* 2 (K). 
We denote by I 2 the space of sequences (ujk)(j,k)&z 2 with the norm \\(ujk)(j^k)&(z 2 )\\ 2 (2 = 
E(jfc)e(z 2 ) 2 2t:? |%'fc| 2 - One caracteristic that plays an important role in operator approxi- 
mation is the semi-orthogonality coefficients B t ,b t > such that y(ujk)j,k&z £ ^t> 

ht/l II /jUjkipjkWH* < II "ifc^'fell-ff* - 5 * /~! II / I UjklpjkWH* 



Shannon wavelets have this particularity to have perfect low-pass and high-pass 
filters * 

1 if f € [-| +2/C7T, § + 2fe7r], fc G Z 



m(f) 
n(0 = 



if £e [f + 2kn, ^ + 2/cvr], fce 



2 i , 2 

-e - * if £e [f + 2fcvr, ^ + 2fcvr], 
if £e [-| + 2/C7T, f +2kn], k£Z 

Then the corresponding scaling function writes: 



1 



m(0 



n(0 



m(0 m(g) 
»(€) 



0(0 



V>(2£) 



2n 



Figure 1: Representation of the compact support of the Shannon filters and of the Fourier 
transform of the wavelets. The vertical axis has no particular meaning, but the bars 
represent the compact support of the functions. 

0(0 = X[-7r,7r](£) , = 

TTX 

and the wavelet: 

~ _i£/2 m ; / \ sin27r(a; - 1/2) sin7r(:r - 1/2) 
m = e XH^ M ,r,*r](0 , = 27r(x _ 1/2) ^ _ 1/2) 

where % stands for the characteristic function i.e. xe( x ) = 1 if x £ E,0 if x tfi E. 

In the multidimensional case, the tensorial Shannon decomposition can be written 
as follows: 

Let u : M. d — > R m . The Shannon decomposition of u is given by: 

u=E u J (3.1) 
jez d 

with 

d 

supp uj C JJ[-2 J ' i+1 vr, -2 ji ir] U [2 ji ir, 2 jl+1 vr] 
i=i 

For each scale parameter j € Z d , and for each component I = 1 . . . m, we have: 
n,j(x) = £ 2j/ 2 ^j k Vf(2 J1 x! - . . . ^(2*a* - h) . . . ifi(2?*x d - k d ) 

where |j| = Y2i3i an d ipf are wavelets of Shannon type, i.e. supp^ C [— 2ir, —ir] U [ir, 2ir]. 
3.1 Shannon wavelet packets 

With the above filters m(£) and n(£), we can define the Shannon wavelet packets. The 
wavelet packets are defined by applying the filters m and n to the wavelets. Hence we 
obtain two new wavelets tpm) and V'(io) that are twice better localised in the Fourier space 
(i.e. the compact supports of their Fourier transforms are twice smaller): 

^(2£) = m(0m (3.2) 
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^(2£) = n(OHO (3.3) 

The two of them are necessary to expend L 2 (IR), i.e. L 2 (M) = span{tpnQ\^ (2^x — 
2k), V'(ii) (2- J ic — 2k)}j t k£Z- More precisely, the wavelet space at level j, Wj, admits 
{V'(io)) (2 J x — 2k), V'(ii) ~ 2&)}fcez as a Riesz basis. 



*K0 



^(25) 



m(£) m(,;) n(£) m(£) 



i i 

2tt 



$(10) (0 ^(11) (?) 



Wavelet 
packets 



Figure 2: The construction of Shannon wavelet packets represented in the Fourier space 

The operations (3.2) and (3.3) on the wavelets can be iterated as many times as one 
wants, and the Fourier support can be shrunk as desired. In practice this operation can 
also be applied to usual wavelets but doesn't come out with good results. Getting a better 
frequency localisation for usual wavelet packets is still a challenging problem. 

The Shannon wavelet packet decomposition gives us more oportunities to approximate 
any operators than the classical Shannon decomposition (3.1). We can report to part 8 
for their actual use. 

4 Solution of elliptic PDE's with wavelets — the Richardson iteration 

In their paper [3], A. Cohen, W. Dahmen and R. DeVore consider a simple method to 
solve elliptic operator equations with wavelets. Let first consider m = n = 1. In order to 
find the solution of the differential equation: 

Au = v (4.1) 

where A is a linear differential operator and u the unknown function, they use the ex- 
pansions of u and v in wavelet bases. We denote by u = {ujk)jk the vector of wavelet 
coefficients: Ujk =< u,ipjk > with {ipjk}j,k and {tpjk}j,k two dual wavelet bases. Then the 
expansion of u writes u = Y2j,k Ujktpjk- 

Let A be the variational discretisation of A expressed in the wavelet basis {ipjk}(j,k) 
(it is called the Petrov-Galerkin stiffness matrix) A = (< A^jk,ipj>k' >)j,j',k,k' and D a 
wavelet preconditionner associated to the wavelet expansions (usually, it is the diagonal 
of A, that has the form Diag(2* J )). We assume that A is continuous from H 1 / 2 to i7~*/ 2 
and a =< A-, - > is coercive. 

In order to diversify the considered wavelet transforms in dimension d > 1, let J denote 
the scale indice set (J countable, e.g. Z x {0, l} d * for the MRA where * means deprivated 
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of the element (0, . . . , 0), and J = % d for tensorial wavelets). For a function a : J — > R we 
define the £ 2 norm on the wavelet coefficients u = (wjk)j e j 5 kez d by 

ll(«jkW )kez -lk= E 2Q(j Vj k | 2 

jG J, k€Z d 

Usually, the t\ corresponds to either the case a(j) = tji in the MRA case j = (ji,e) G 
Z x {0, l} d * , either a(j) = f max(ji, . . . for the tensorial wavelets. 

Then L> _1 is continuous from £ 2 _ t i 2 to (see [12]). We write the sequence (u n ), 
thanks to a Richardson iteration associated with a multiscale preconditionnning, starting 
from u Q = : 

=Mn-l (4-2) 

Then this method is said to converge if 

3/9 < 1, \\V -AuJp t/2 < p\\v - Mn-l\\p_ t/2 

As we have: 

v ~ Au n = (Id- AD~ 1 )(v - Au^i) 
the algorithm converges if p = \\Id — ^^4^0 1 1 1 < 1, in the operator norm. And 

U n —^n-too u with A U = V 

That is, as u = Yljk u jk'4'jk G if*/ 2 , 

Au = v 



From now on we think of A as being an operator with constant coefficients. Let m, n 
be unspecified natural numbers. Hence we switch to vector spaces. If we denote by M(£) 
the symbol associated to A, we can express A after a Fourier transform of the equation 
(4.1) as 

M(£)u = v and the pseudo-inverse solution u = M(£)^v 

with u G (i7'/ 2 (M d )) m , v G (H^ 2 (R d )) n , Af(f) G C nxm and M(£) f the pseudo-inverse of 

M(e). 

Remark that if m = n and M(£)M(£)t = Id, M(£) f = Af (f) -1 . 

The idea for solving = v is the following: we decompose v in a wavelet basis that 
splits the support of v 

V = E V J 

If we denote by (ei, . . . ,e n ) the canonical basis of W 1 then vj is the projection of v in 
the wavelet level W- } = span({\&ijk}kez d > • • • ' i^njk}k£Z d ) with each component vgee of 
v decomposed in the wavelet basis jk}jej,kgz d (further we'll need this generalisation 
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of [3] which uses an MRA). 

For example in the tensorial case, we have J = 7L d and: 



w/>iifci(zi) • • • i>j d k d ( x d) 



This modification of the MRA case will prove usefull in part 10. 

Let us assume that for each j G J, vj and uj are compactly supported (wavelet decompo- 
sitions give us the opportunity to do this with the desired accuracy). For each j € J, we 
build a matrix M u . (M^. G R nxm depending on the compact support of vj) such that: 



M(£) for £ G supp(vj) 



(4.3) 



Then we approximate the relation M(£)u = v by: 



Vj eJ, Vke Z d , M u 







U ljk 









(4.4) 



In the view of Richardson iteration, we take as a preconditionnner D = M^.Pj, 
where Pj is the projector on the wavelet level Wy Then the corresponding discrete pre- 
conditionner D which applies to wavelet coefficients is D = X^jeJ ^^j—i wnere LZj is a 
diagonal matrix with ones on the lines (k, j) keZ d and zeros everywhere else. In the case of 
tensorial wavelet basis (J = Z rf ), the space Wj is the L 2 closure of the space generated by 
the family 

{*Ajk}i<^< m ,keZ d = {WVifciX- • ,x ^idfcj2irii2^ k6ZdU " ' • u {(0 1 —_ 1 0,V'jifciX- • -x^JIkez" 

m— 1 m— 1 

In the following, we use the notation = ^ M Wj i^. 

If we write the sequence (4.2) with v n = v — Au n , it comes: 

u = , v = v , u n+1 = u n + M+V n and v n+1 = v n - A(u n+1 - uj (4.5) 

Theorem 4.1 Let M(£) 6e t/ie symbol matrix associated to A : hi 1 ! 2 — > H^ 1 ! 2 continuous. 
If the wavelet basis {*4jk}i<£< T ra,jeJ,keZ< i provides a Riesz basis of H ±l / 2 (i.e. the asso- 
ciated decompositions vhv, h ±1 / 2 _ > and reconstructions vh v, ^±*/2 — * H^ 1 / 2 
are continuous). 

Moreover, we suppose we have constructed for all j G J matrices M w . G R nxm s-uc/t £/iaf 



^•M(i.i^ : iJ */ 2 — > ff*/ 2 is continuous. We also assume that the wavelet decom- 



position v i— > (vj)j 6 j satisfies: 

3B>0 such that Vv G #~ </2 , ||(Id- AMj, )v||^_ t/2 < ||(Id- AMj, )vj||^. 

jeJ 

7/ i/iere exisi a rea/ number p > suc/i £/iat: 

VjGJ, |||(7d-A^W|||<p 



fc/2 
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i.e. 

Vj G J, Vvj G Wj, IKJd-AA^JVjl^-t/a < p||yj|| H _ t/a 

i/ten /or p smaH enough, the sequence (u n ) n eN defined by (4-5) converges in £ 2 , 2 to the 
wavelet expansion u of u s-uc/t i/iai: 

A u = v and we have Au = v 

proof: 

The graph of continuous operators can be summarized as follows: 

wavelet transform 
u G H 1 / 2 * ► u G i 2 t/2 

A[ ]Ml Ml] I A 

v = A u G H~ 1 / 2 < > v = iu G 

wavelet transform 

The operator is not the inverse of A but its approximation. 
As the wavelet decompositions are continuous, 

3b,B>0 such that Vver'/ 2 , 6^ ||vj||^_ t/2 < ||v||^_ t/2 < £^ ||vj||^_ t/a 

jeJ jeJ 
When 6 = B = 1, the wavelet basis is said to be semi-orthogonal. 
Then we have: 



,B , 



|2 

'n\\ff— 1/2 



jeJ jeJ 

If p 2 B/b < 1, as Mj) is continuous, the serie ^ra-^i v « converges in the Banach space 
H 1 / 2 to a solution u of the equation Au = v. 

The ideal wavelets that provide a minimal compact support for the Fourier transform 
are the Shannon wavelets. In this the compact supports of the Fourier transforms 

of wavelets from different levels are disjoint, we have B = b = 1 for all A. Shannon 
wavelets have an infinite support and are not used in practice. But, in first approximation, 
all wavelets behave as Shannon wavelets with more or less accuracy. 

Remark 4.1 In the case of Shannon wavelets, as 

9 |U f supp(*i) = E 

the equation 4-4 * s equivalent to 

VjGJ, M wj u(0=v(f) for ee^supp(*i) 

In the futur, that will allow us to express this relation using the components vj of the 
Shannon decomposition of v. 

One can also remark that as M u doesn't depend on £, we can apply the operator Mjj in 
the physical space (expressed with wavelets) and not in the Fourier space. 
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5 Multiresolution analysis (MRA) versus tensorial basis 

There are two main different kinds of wavelet decompositions for a function on M. d with 
d>2. It can be decomposed either in a multidimensional multiresolution analysis or in a 
tensorial basis. In an MRA, the wavelet decomposition of a function u in 2D writes: 



u{xi,x 2 ) = X! Yl u f,kiM ^Oj.fcxO^i) Vij,k 2 ( x 2) + Y u f,kiM Wi.fciOci) ^li.*a( 
jez \(fc 1 ,fc 2 )ez 2 (£.-i,fc 2 )ez 2 

(fci,fca)eZ 2 

where we used the notation ipj^(x) = ^{2? x — k). 
While, in a tensorial basis it writes: 

u(x 1 ,x 2 ) = YY Y u h,hMM V ; o(2 il a;i - h) i/j 1 (2 j2 x 2 - fc 2 ) 
iieZi 2 eZ(fc 1 ,fc 2 )ez 2 

These two decompositions correspond to two different partitions of the Fourier space (i.e. 
frequency domain). Both of them are represented in figures 3 and 4. On each figure, 
in the last square, which corresponds to the wavelet transform, the low frequencies are 
localised in the upper left corner of the square of coefficients, and the high frequencies in 
the bottom right. 



$J-l,k 




^(1,0) 























< 



High frequencies 



Figure 3: Splitting of the Fourier modes induced by the 2D-MRA wavelet decomposition 



5.1 Convergence theorems with Shannon wavelets 

To begin with, we'll only consider approximation matrices that are constant over each 
frequency domain indexed by j GZf 

The two previous decompositions induce different conditions for the approximation of 
the matrix M(£). Adding to part 4, we distinguish MRA and tensorial wavelet convergence 
theorems as follows: 

Theorem 5.1 (MRA) If the symbol matix M(£) admits a pseudo-inverse M(£)t such 
that M(f )M(f ) f = Id V£ ^ (0, . . . , 0), and if 3p < 1 such that Vj G Z and Ve G {0, l} d \ 
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<pJM x i'JM 



Low frequencies 



x 4>j,k 2 (y) 









































).Cfcl.*2j 



High frequencies 



Figure 4: Fourier splitting induced by the tensorial wavelet decomposition 



{(0,...,0)} 3M u , a e 



such thatV£ E nf=i ±[ei2%, (^ + 1)2%], \\Id-M(£)M t 



< 

p then the sequence (4-2) using the MRA decomposition with Shannon wavelets, converges, 
proof: 

We recall that we are in the case b = B = 1 of theorem 4.1 since we deal with Shannon 
wavelets. 

The partition of the support of v operated by the MRA decomposition is the following: 

J = {(j, £ )eZx{o,i}^} 



V = V J With SU PP 9 i C ±[ £ i 2i7r > ( e i + i) 2 ^] 

This writting is due to the fact that the case £j = corresponds to a scaling function (fij 
for the variable Xj, and £j = 1 to a wavelet function ipj. Owing the fact that supp^- C 
[— 2^ ir, 2 J 7r] and suppV'j C [— 2- ?+1 7r, — 2 J 7r] U [2 3 V, 2- j+1 7t], we obtain the sets indicated in 
the theorem. This case is represented in figure 3. 
Then we apply theorem 4.1 to obtain the convergence. 



Remark 5.1 The fact that M(£) admits a pseudo-inverse M(£)t such that M(£)M(£) 
Id is implied by 3p < 1, 3M W . £ M nxm suc/« i/tai \\Id - M(£)Ml \\ < p. 



Theorem 5.2 (Tensorial wavelets) If the symbol matix M(£) admits a pseudo-inverse 
M(£)t s-uc/j that M(£)M0 = Id V£ E M d \ {£ E M d sucfe tto* &...& = 0}, and 3p < 1 
smc/i that Vj E Z d 3Af Wj € R" xm suc/i V£ E nf=i ±[2%, 2* +1 7r], M^M^ || < 

p then the method converges with the tensorial wavelet decomposition. 



proof: 

Anew, we use theorem 4.1 with J 



and 



v = with suppVj C ±[2 ji ir, 2 ji+1 ir] 

jeJ 



This wavelet decomposition part the frequency domain as represented in figure 4. 
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Remark 5.2 If we consider only constant matrices operating on the wavelet coefficients, 
the resulting operations on the Fourier transform of the functions are symetric by reflection 
along all axes: 

Vie{l,...,d}, Mufo,..., -&,...,&) =M W (&,... ,&,...,&) (5.1) 
On the other hand, as we deal with real functions, the approximation matrix must be real. 
Remark 5.3 The best aproximation M^. of M(£), for the inversion is given by 

M u . = arg min sup 

Example 5.1 As we'll see in part 9, the operator A -1 matches the two cases. Wavelet 
algorithms converge in the MRA context and in the tensorial one. 

Example 5.2 If we consider the 1-D symbols p : R — > C, £ i— ► that are continuous, an 
example that doesn't match the conditions of the theorem is given by a symbol p such that: 
p(0 = e ^ f or £ £ I 71 "' 27r]. TTien i/ie operator whose symbol is p can't satisfy 3p < 1, £ R 
swc/i thatV£ £ [7T, 27r], ||1 — p(£)^~ 1 || < P- 

Example 5.3 On i/ie other hand, the 1-D symbols p : R — > R, £ h- ► that are contin- 
uous, verify V£ 7^ 0, p(£) 7^ and 

SUp €e r2»V l2 J+%](|p(0l) 

sup w Pi /fin < +°° 

can 6e approximated by a constant u>j on each interval ±[2-%, 2- J+1 7r] mi/i optimal value 
verifying 

uj 1 = \( m 7 1 + M j 1 ) 

with mj = inf ee[2 i 7r)2 ,+i 7r] (|p(£)|) and M,- = sup fe[ 2j ff)2 i+i w ](|p(0|) 

77ia£ is the case for functions with polynomial increase, since for £ a , ^j^ja = 2 Q . 

Example 5.4 In 2-D, even for real operator matrices, the approximation by constant 
matrices can fail. For instance, if we consider a symbol matrix M such that: 

any wavelet approximation by constant matrices \x £ R 2x2 fails: either \\Id + fj,\\ > 1 or 
H-fd-HI > 1- 

6 Derivation of wavelets 

P. G. Lemarie-Rieusset [14] showed that derivating or integrating a biorthogonal wavelet 
basis provided a new wavelet basis. It allows us to construct two different one-dimensional 
multiresolution analyses of L 2 (R) related by differentiation and integration. 
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Theorem 6.1 (Derivation of wavelets) [14J Let (V^)j & z be a one- dimensional MRA, 
with a differentiable scaling function ip\, (Vq 1 = span{(pi(x — k), i: £ Z}J, and a wavelet 
There exists a second MRA (V^ )jez with a scaling function (fo (V® = span{(fo(x — k),k € 
Z} ) and a wavelet tpo satisfying: 

ip[(x) = ip (x) -ip (x- 1) il>[(x) = 4 iI>q(x) (6.1) 

Expressed with its Fourier transform this relation writes: 

mMO = 4^(0 

The filters (mo, m^) and (mi, m*) attached respectively to the MRA 's (Vf)jeZ an d (^j 1 )jez 
verify: 

2 1 + e*£ 

m o(0 = ~. : — —7 m i(0 and m o(0 = — o — m i(^) 
1 + e l? z 

If the wavelet ipi is C n and has p zero moments (i.e. ip\ is derivable p — 1 times in a 
neighborhood of and '0^ (0) = for < k < p — 1) after such an operation, the wavelet 
tpo has regularity C n_1 and p + 1 zero moments. 

Remark 6.1 As the Shannon wavelets are C°° and have an infinite number of zero mo- 
ments, they can be derivated or integrated in order to obtain biorthogonal wavelets sat- 
isfying the relations (6.1) of the theorem 6.1. And we can iterate the derivation or the 
integration of these wavelets as many times as we wish in order to obtain derivatives of 
arbitrary order: . . . , tp-2, tp-i, ■ ■ ■ , 1P2, ■ ■ ■ , with ipo the original Shannon wavelet. 

On account of the above remark, we can introduce a new operation thanks to the wavelet 
decomposition of a function v. Indeed, if we use the tensorial wavelet decomposition, 
we can derivate or integrate in every directions. For instance, if we write the wavelet 
decomposition of v with wavelets V'o f° r each tensorial components except for i for which 
we take ipi where ipo and ip\ are related by the derivation relation (6.1) as in theorem 6.1. 

v(x)= d j ^o^ il ^i-ki)-..M^^-k i ).--M^ jd Xd-k d ) 

Then if we put for u: 

u(x)= ^ 4-2-^ jk Vo(2 il x 1 -k l )...^{2 jl x i -k i )...^ d x d - k d ) 

We obtain: 

Qv ^ ^ 

u(x) = — — (x) or, in Fourier u(^) = iCiv(0 

OXi 

7 Constructible approximations 

Here we restict ourselves to the case when m = n. From the results of the previous section, 
it comes: 
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Theorem 7.1 (Set of constructible operators) The set of symbol matrices that are 
constructible by multiplying the wavelet coefficients by some constants depending on the 
parameter j and by derivating wavelets as in theorem 6.1 is the M-algebra of C nxn (£) 
generated by the elements {(5i,j)}i<i,j<n, {i£il}i<i<d an d I}i<i<d) with (Sij) denot- 
ing the matrix which is zero everywhere except at line i and column j where it is 1. 

This theorem enables us to diversify our wavelet approximations of differential operators 
and extends the result of section 4. But it still remains rather limited since for instance, 
in dimension larger than 2, we cannot reach A -1 , the inverse Laplacian, with these oper- 
ations. 

8 Convergence with the wavelet packets 

Let A be an operator from (H t / 2 (R d )) n to (H~ t / 2 (R d )) n , having a continuous symbol 
M(£) almost everywhere invertible on W 1 in the sens of Riemann measure (i.e. such that 
for all compact sets K of R d , K n (det(M))" 1 ({0}) the subset of K where M(£) is not 
invertible has a vanishing Riemann measure), and verifies the condition (5.1). Then it can 
be approximated by constant matrices M u . with Shannon wavelet packets providing an 
ad hoc partition of the frequency domain. 

Theorem 8.1 For all linear operator A satisfying the above conditions, we can numeri- 
cally solve the equation Au = v with a wavelet packet method, i.e. Ve > 0, we can find u e 
such that || v — Au t \\ < e, thanks to the wavelet algorithms described in part 4- 

Proof: First we build a finite set of rectangles {cjj}j g j of the type 

d 

u; j = n [2^,2^ + 1)] 

i=l 

withj = (ji, . . .,j d Ji, ■ ■ ■ ,&d), ji G z and h G N, such that if ft = Ujgj w j tlien < 
e/2, and for $ = (2^4, • • • , 2*4, . . . , 2?H d ) G wj, 

sup [ sup 1 1 Id — M(£)M -1 (£j)|| j < 1 
J J 

Remark that, as a consequence, M(£) is invertible on 0. 

On the set 0, we can apply the wavelet algotithm of part 4 with Shannon wavelet packets, 
and get u e such that \\v\q — A\qu £ iq\\ < e/2, then we extend u e to M d by taking u e \cQ = 0. 

Remark 8.1 This approach is valid for Shannon wavelets. But in practice we would like 
to use other wavelets for which the frequency partition induced by the wavelet packets is 
very difficult to control. 

Remark 8.2 Nevertheless, this method should help us to improve the existing algorithms 
(see parts 9 and 10). 
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9 Implicit Laplacian 

First, let us begin with the elliptic operator A = (Id — a A), a > 0. The above study allows 
us to state precise properties for the wavelet iterative algorithm solving (Id— aA)u = v 
(in Navier-Stokes a = v5t). 

First, we consider a Shannon wavelet decomposition for u and v: 

(id - «a) u j = E v j 



then, for each j G Z rf , we solve 

(Id - aA) Uj = Vj 
with supp(uj) = supp(vj) C nf=i ±[2^ +1 7r,2^7r] 

We approximate M(£) = (l+a\t\ 2 ) Id with £ G flti ±[2 il+1 vr, 2*tt] by M Wj = (1+au?) Id 
with wj £ 1 properly chosen. 

Now, we'll see what is the most convenient value for the parameter ujy The symbol 
matrix (Id - M(£)M~ l ) of the operator (Id - AM' 1 ) is a diagonal matrix Xld whose 
element A 6 1 is equal to: 

1 + auj? 

a?+b? 

This expression has a minimal maximum for |£| 2 G [a? , 6?] at = J 2 J ■ 
If we take this value for Wj, then 

Ve G M d / o? < Kl 2 < &?, |A(/d - m^m- 1 )! < 



2 + a(6? + a?) 

Then in the tensorial case, with the Shannon wavelets, as for fixed j, Vi G {1, . . . , d}, G 
±[2 J 7r, 2 J+1 7r], then 6j = 2a- } and |A| < 2a r^+5a ' ^' ie wors ^ case a PP ears f° r a ~^ +oo, 
then the convergence rate tends to p = I = 0.6. 

In the case of the MRA, we have always b?- < (d + l) a ijj that allows the algorithm to 

converge too: p = 2g -^ d+2)a ^a- +0 o 3^. 

If we now consider what should be obtained with the Shannon wavelet packets, we 
take b = 3/2 a in the tensorial case, then if a — > 00, p = ^ (~ 0.38). Roughly speaking, 
we improve the convergence by a factor 2 for each wavelet packet refinement. 

Observed convergence for the implicit Laplacian 

For the operator Id — a A, the convergence of the wavelet algorithm is very fast if a is 
small compared to the smallest computed scale. Anyway, for the Laplacian operator A 
(a — > +00) the same algorithm still converges, and the observed convergence rate is around 
0.75 for spline wavelets of order 4. 



15 



10 Leray projector 

Principle of the Helmholtz decomposition: 

Let u S (L 2 (W 1 )) n be a vector field. We can decompose u as follows: 

u = u div + u cur i where u div = curl tp , u curl = Vp 
The functions curl tp and Vp are orthogonal in (L 2 (M. n )) n and are unique. 

(L 2 (m n )) n = H div0 (M") X H curl0 (R re ) 
In Navier-Stokes, this decomposition is very important to project the term u.Vu onto 

H 

divo(^ n ) the space of divergence free vector functions (see [1]). 

The wavelet algorithm that will be studied here was originally designed in [8]. The 
proof of its convergence for the Shannon wavelet is new. The wavelet iterative al- 
gorithm for P, the orthogonal projector (in L 2 ) on H d i v o(P"), uses the result on the 
wavelet derivation of part 6. We approximate the Leray projector P: Pu = M(£)u with 



M(f ) = Id 



1 

W 



x [£i...£„] 



by 



Id 



\ 



x &...£„] 



Id 



\u\ 



J 



6 



where we used the notation uj~ = 2 3i . 



Remark 10.1 One can notice that this projector M u actually projects the vector field u 
on the space of divergence free vector functions (indeed < [£i . . . £ n ], M u u(^) >= 0). 

We also have to extract the gradient part of u by the approximation of VA _1 ( div-): 

£l 



I, 



1 



uj 



Theorem 10.1 The sequence (4-5) writes: 

v = u and v n+ i 



£i " ' & 



flO.ll 



27ms sequence (v n ) goes to m -L 2 norm in the case of Shannon wavelets. The iteration 
of the operation (10.1) allows us to separate the divergence free part u div = ^2 ne ^M u v n 
of the function u from its gradient part u CUI \ = J2neN v n . 
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256 x 256 points 
2 zero moments 

512 x 512 points 
2 zero moments 



1024 x 1024 point 
2 zero moments 



256 x 256 points 
3 zero moments 



Number of iterations 



Figure 5: Observed convergence rates for various wavelets and numbers of grid points 



proof: 

The matrix (Id - M u - L w ) can be written, for £ G [Ji ±[2 jl vr, 2?< +1 7r] 



(Id-Mu-Lc) 



6 



X [6 . . . Cn] X (Id - L u ) 



All its eigenvalues are equal to zero except one that equals 

A ^ 



X = 1 



i.e., with Cfc = ^, 



n o 



\k=l 



Then the Kantorovitch inequality yields 



IAI < 



1 / min fc |Cfc| . max fc |Cfc 



+ 



1 



4 \max k |Cfc| min fc \( k \ 
If |Cfc| S [afc, 6fe] for each k, we put a = min^ a& and 6 = max^ b k . Then 



1 (a b 
|A| < - ( - + - 



4 V b 



1 



For the Shannon wavelets, b = 2a, then convergence is assured since |A| < ^ (~ 0.56). 
In the case of the Shannon wavelet packets, b = |a, the convergence rate should be 



|A|< 



25 
144 



0.17). 



Observed convergence for the Leray projector 

The convergence has been tested successfully on variate 2D and 3D fields. The observed 
convergence rates with spline wavelets of order 2 and 3 are around 0.5 (see figure 5). The 
reference [8] provides technical explanations for the impementation of this algorithm. 
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11 Conclusion, perspectives 

This work provides an original point of view on the wavelet algorithms and links two 
achievements of the wavelet theory: 

• the wavelet approximation of operators like it is presented in A. Cohen's, W. Dah- 
men's and R. DeVore's work [2, 3, 4, 6], 

• the divergence-free wavelet transform derived from P.-G. Lemarie's and K. Urban's 
works [14, 17] as well as the author's [7, 8, 9]. 

The main results of this paper are the establishment of general conditions for the conver- 
gence of wavelet algorithms with Shannon wavelets, the theoretical construction of wavelet 
algorithms in order to approximate operators with constant coefficients, the exact compu- 
tation of the convergence rates and their optimisation for the implicit Laplacian operator 
and the Leray projector. It also gives a simple view of the wavelet preconditionning and 
a glance to what could be done thanks to wavelet packets. Indeed, wavelet packets seem 
to provide a powerfull solver for many kinds of PDE's but their use is still theoretical and 
prospective. 

The progress to achieve should go toward a more general frame for the proofs of the 
convergence of wavelet algorithms that would allow us not to restrict ourselves to Shannon 
wavelets. As a result we should be able to introduce wavelets on the interval in this frame. 
We are still missing efficient wavelet packets concerning the frequency localisation. An 
interesting perspective would be to study what happens in the case of operators with non 
constant coefficients 

These theoretical assertions have already an application since these technics are de- 
rived to simulate the Navier-Stokes equations (see [7], and a forthcoming paper in Siam 
Multiscale Simulation and Analysis). 
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